if !d.name eq 'PS' then begin
    device,xsize=9,ysize=15,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end

@eu_read_bcprofl
;!p.multi=[0,3,2]
;!p.charsize=2.5
!p.multi=[0,1,1] & !p.charsize=1.0 
!x.margin=[7,0] & !y.margin=[2,3]

default,nn,5
ntimes = (size(u))[3]
good = (ntimes-nn)+indgen(nn)-1
print,'doing averages over',nn, ' xaver records ...'

ruu = (total((total(u2,1)/m),1)/l)-(total((total(u^2,1)/m),1)/l)
rvv = (total((total(v2,1)/m),1)/l)-(total((total(v^2,1)/m),1)/l)
rww = (total((total(w2,1)/m),1)/l)-(total((total(w^2,1)/m),1)/l)


uu=fltarr(l,m)
vv=fltarr(l,m)
ww=fltarr(l,m)
oox=fltarr(l,m)
ooy=fltarr(l,m)
ooz=fltarr(l,m)
pt0=fltarr(l,m)   ;potential temperature
pre=fltarr(l,m)   ;pressure

omega=fltarr(l,m)
ux=fltarr(l,m)
uy=fltarr(l,m)

quu=fltarr(l,m)
qvv=fltarr(l,m)
qww=fltarr(l,m)
qwv=fltarr(l,m)
qwu=fltarr(l,m)
qvu=fltarr(l,m)
rsinth=fltarr(l,m)
rurms2=fltarr(l,m)
frmc=fltarr(l,m)
frrs=fltarr(l,m)
fthmc=fltarr(l,m)
fthrs=fltarr(l,m)

frt=fltarr(l,m)
ftht=fltarr(l,m)
fx=fltarr(l,m)
fy=fltarr(l,m)

urms2t = ruu+rvv+rww
urms2 = (total(urms2t[good]/float((size(good))[1])))

; reynolds stresses
for k=0, l-1 do begin
  for j=0,m-1 do begin
    rsinth[k,j]=rr*r[k]*sin(theta[j])
    uu[k,j] = total(u[j,k,good])/float((size(good))[1])
    vv[k,j] = total(v[j,k,good])/float((size(good))[1])
    ww[k,j] = total(w[j,k,good])/float((size(good))[1])
    oox[k,j] = total(ox[j,k,good])/float((size(good))[1])
    ooy[k,j] = total(oy[j,k,good])/float((size(good))[1])
    ooz[k,j] = total(oz[j,k,good])/float((size(good))[1])
    pt0[k,j] = total(ptemp[j,k,good])/float((size(good))[1])
    pre[k,j] = total(p[j,k,good])/float((size(good))[1])
    quu[k,j]=rho[j,k]*total(u2[j,k,good]- u[j,k,good]^2)/float((size(good))[1])
    qvv[k,j]=rho[j,k]*total(v2[j,k,good]- v[j,k,good]^2)/float((size(good))[1])
    qww[k,j]=rho[j,k]*total(w2[j,k,good]- w[j,k,good]^2)/float((size(good))[1])
    qwv[k,j]=total(rwv[j,k,good]- rho[j,k]*oz[j,k,good]*v[j,k,good])/float((size(good))[1])
    qwu[k,j]=total(rwu[j,k,good]- rho[j,k]*oz[j,k,good]*u[j,k,good])/float((size(good))[1])
    qvu[k,j]=total(rvu[j,k,good]- rho[j,k]*oy[j,k,good]*u[j,k,good])/float((size(good))[1])
    rurms2[k,j]=quu[k,j]+qvv[k,j]+qww[k,j]
    ux[k,j]=ww[k,j]*sin(theta[j])+vv[k,j]*cos(theta[j])
    uy[k,j]=ww[k,j]*cos(theta[j])-vv[k,j]*sin(theta[j])
; r components of the angular momentum flux. all normalized with urms2
    frmc[k,j]=rsinth[k,j]*rho[j,k]*uu[k,j]*ww[k,j]/rurms2[k,j]
    frrs[k,j]=rsinth[k,j]*qwu[k,j]/rurms2[k,j]
; th components of the angular momentum flux. all normalized with urms2
    fthmc[k,j] = rsinth[k,j]*rho[j,k]*(uu[k,j])*vv[k,j]/rurms2[k,j]
    fthrs[k,j] = rsinth[k,j]*qwv[k,j]/rurms2[k,j]
    frt[k,j] = (frmc[k,j]+frrs[k,j])
    ftht[k,j] = (fthmc[k,j]+fthrs[k,j])

    fx[k,j]=frt[k,j]*sin(theta[j]) + ftht[k,j]*cos(theta[j])
    fy[k,j]=frt[k,j]*cos(theta[j]) - ftht[k,j]*sin(theta[j])

  endfor
endfor

omega_0=1.e9/(28.*24.*3600.)
omega[*,1:m-2]=1e9*(uu[*,1:m-2]/(2.*!pi*rsinth[*,1:m-2]))
omega=omega+omega_0
print,"Omega_0=",omega_0
nlev=64
thg=where(th GE -89 AND th LT 89)

pnorm=1.71889e+13
thnorm=3.4976e6

;print,'levels Omega = ',minmax(omega)
;lev = 2.*max(abs(pre/pnorm))*indgen(nlev)/float(nlev-1)-max(abs(pre/pnorm))
;contour,pre[*,thg]/pnorm,x[*,thg],y[*,thg],/fi,nl=64,/iso,$
;	title='!8p`/p!D!6ref!N',lev=lev,xtitle='!8x!6',ytitle='!8y!6'
;colorbar,range=[min(lev),max(lev)],pos=[0.22,0.4,0.24,0.68],/vert,ytickformat='(E10.1)',yticks=2,ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],yaxis=0,char=1.1
;;partvelvec,ux[*,thg],uy[*,thg],x[*,thg],y[*,thg],/over,fraction=0.08,length=0.08
;plots,circle(0,0,0.98),thick=3
;plots,circle(0,0,0.62),thick=3
;plots,circle(0,0,0.71),thick=2,li=2
;xyouts,0.8,0.8,'!6A.',charsize=2

;
;for k=0, l-1 do begin
;  for j=0,m-1 do begin
;    if (pt0[k,j] le 0) then begin 
;      pt0[k,j] = -alog(abs(pt0[k,j]))
;    endif else begin
;      pt0[k,j] = -alog(pt0[k,j])
;    endelse
;  endfor
;endfor

lev = 2.*max(abs(pt0/thnorm))*indgen(nlev)/float(nlev-1)-max(abs(pt0/thnorm))
lev = 2.*max(abs(pt0))*indgen(nlev)/float(nlev-1)-max(abs(pt0))
print,'Total radial flux = ',minmax(lev)
contour,pt0[*,thg],x[*,thg],y[*,thg],/fi,nl=64,/iso,levels=lev,$
	title='<!7h`>!6!N',xtitle='!8x!6',xstyle=4,ystyle=4

colorbar,range=[min(lev),max(lev)],pos=[0.27,0.35,0.295,0.65],/vert, $
	ytickformat='(E10.1)',yticks=2,ytickv=[min(lev),0.,max(lev)],yaxis=0,char=1.1
plots,circle(0,0,0.96),thick=3
plots,circle(0,0,0.602),thick=3
plots,circle(0,0,0.718),thick=2,li=2
oplot,[0,0],[-0.96,-0.602],thick=3,li=0
oplot,[0,0],[0.602,0.96],thick=3,li=0
save,filename="ptemp.sav",pt0

th_mean=fltarr(l,m)
th_ave=total(pt0,2)/double(m)
for k=0, l-1 do begin
  for j=0,m-1 do begin
   th_mean[k,j]=pt0[k,j]-th_ave[k]
  endfor
endfor

print,';;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;'
;
!p.multi=0
end
